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I.  INTRODUCTION 


ISEAS  (Interactive  Simulation  of  Engagements  At  Sea)  is  a  computer 
program  currently  under  developement  at  the  Naval  Postgraduate  School 
to  simulate  interactions  between  a  ship  loaded  with  surface-to-air 
missile  systems  and  an  attacking  aircraft  or  anti-ship  missile.  The 
portion  of  ISEAS  described  in  this  thesis  is  a  translation  of  the 
Fortran  program  MICE  II  [Ref.  1 : pp  4-20]  for  the  flyout  of  a  surface 
to  air  missile  to  the  C  language.  Three  different  missile  systems  are 
simulated.  The  C  program  includes  coordinate  transformations,  launcher 
bearing  and  elevation  calculations,  the  proportional  navigation  equa¬ 
tion  for  missile  guidance,  and  the  missile’s  equations  of  motion.  The 
program  contains  a  description  of  the  missile  engagements  with  an  air 
target  and  provides  a  basis  for  the  evaluation  of  surf ace-to-air 
missile  launching  opportunities  and  subsequent  missile  performance 
under  realistic  encounter  conditions.  This  missile  program  provides 
the  basic  framework  from  which  further  simulations  of  increased  comp¬ 
lexity  and  sophistication  can  be  easily  implemented. 
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II.  PROPORTIONAL  NAVIGATION 


Guided  missiles  are  navigated  along  a  path  that  leads  to  an  inter¬ 
cept  with  the  air  target.  Proportional  Navigation  is  one  type  of 
navigation  in  which  the  rate  of  change  of  the  missile  heading  is  made 
proportional  to  rate  of  change  of  the  line  of  sight  between  the  missile 
and  the  target.  See  Figure  1  for  the  Proportional  Navigation  Guidance 


Figure  1 .  Proportional  Navigation  Guidance 
Tra jectory 

Trajectory.  The  equation  governing  Proportional  Navigation  is 

6  =  k*0  ( 1  ) 

where 

6  =  time  rate  of  change  of  missile  heading 
(velocity  vector) 

* 

0  =  time  rate  of  change  of  line  of  sight 
k  =  Proportional  Navigation  constant 
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The  missile  rate  of  turn  (0)  is  a  fixed  proportion  of  the  rate  of 
turn  of  the  line  of  sight  (0),  but  the  proportional  constant  itself  (k) 
can  be  variable  during  the  flight.  This  method  is  implemented  in  MICE 
II.  The  input  base  value  k0  ,  is  multiplied  by  the  ratio  of  the  closing 
velocity  to  the  missile  velocity.  Consequently,  due  to  the  changes  in 
velocities,  the  effective  value  of  the  Proportional  Navigation  constant 
is  a  variable  during  the  missile  flight.  There  are  two  special  cases 
of  Proportional  Navigation,  Pursuit  and  Constant  Bearing  or  lead  angle. 
The  Pursuit  Navigation  course  can  be  obtained  from  Proportional  Naviga¬ 
tion  by  setting  k  =  1 ,  and  the  Constant  Bearing  course  can  be  obtained 
from  Proportional  Navigation  by  letting  k  be  very  large. 

A  missile  traveling  a  Proportional  Navigation  course  is  usually 
launched  at  a  lead  angle  which  corresponds  to  the  Constant  Bearing 
path.  (Note  that  a  pure  Pursuit  course  is  actually  the  unique  case  of 
a  Proportional  Navigation  path  with  a  zero  lead  angle  and  a  navigation 
constant  of  one).  Proportional  Navigation  paths  are  less  curved  (fewer 
g’s  on  the  missile)  than  the  Pursuit  course,  but  may  be  more  curved 
than  the  Constant  Bearing  course.  The  advantage  of  the  Proportional 
Navigation  course  is  that  it  provides  a  practical  method  of  causing  an 
intercept  with  a  manuvering  target.  Suitable  values  of  the  navigation 
constant  are  between  three  and  six.  [Ref.  2:p.  67] 

Figure  2  shows  the  two-dimensional  geometry  for  a  Proportional 
Navigation  course.  The  missile  and  target  are  a  distance  r  apart. 
The  line  of  sight  between  the  missile  and  the  target  is  the  angle  0 
relative  to  the  global  reference  (X,Z)  for  all  times  of  flight.  The 
line  of  sight  is  at  an  angle  relative  to  the  missile  velocity  vector, 
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and  the  line  of  sight  is  at  an  angleC^^  relative  to  the  target  velocity 
vector.  The  relative  closing  velocity  between  the  missile  and  the 


target  (r)  can  be  written  as  [Ref.  2: p.  68] 

r  =  V.  *  cos  a:  -  V  *  cos  oC  ^  (2) 

T  t  'rry‘ 

where  *  cos  =  velocity  components  of  the  target  along 

the  line  of  sight 

V  *  cos  0Cm  =  velocity  components  of  the  missile  along 
the  line  of  sight 


The  rotation  rate  of  the  line  of  sight  can  be  written  in  the 
form  [Ref.  2:p.  68] 


r0 


-V, 


sin  cx  +  V  *  sin  oc . 
t  ** 


(3) 


where 


V^*  sin  =  velocity  vector  component  of  the  target 
normal  to  the  line  of  sight 


V  *  sin  Oc^=  velocity  vector  component  of  the  missile 
normal  to  the  line  of  sight 

0  =  rate  of  change  of  the  line  of  sight 

angle 


Equations  (2)  and  (3)  are  the  basic  equations  of  motion  for 
Proportional  Navigation . 
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III.  THE  MISSILE  PROGRAM 


A.  GENERIC  MISSILES 

Simulations  of  three  types  of  surface-to-air  missiles;  Long  Range, 
Medium  Range,  Short  Range,  are  developed  for  ISEAS.  The  general 
default  missile  parameters  for  each  missile  type  are  indicated  below: 


Long  Range 

Medium  Range 

Short  Range 

wt  (kg) 

1066 

581 

231 

body  dia  (m) 

0.35 

0.35 

0.20 

length  (m) 

7.98 

4.47 

3.70 

wing  span  (m) 

0 .91 

0 . 91 

1  .  00 

ref erence 
area  (  sq .  m  ) 

8.78 

4.92 

2 . 33 

range  (km) 

121 

48 

20 

max  alt  (km) 

20 

20 

20 

max  flight 
time  (sec) 

200 

80 

40 

avg  spd  (mach) 

2.5 

2 . 5 

2.5 

thrust  (kN) 

400 

400 

200 

initial 

velocity  (m/sec) 

670 

670 

670 

missile  burn 
rate  (  kg/sec  ) 

3.54 

4.93 

4.67 

drag  coeff 

.2 

.2 

.2 

The  user  will  have  the  option  of  changing  any  of  the  missile  para¬ 
meters  listed  above,  but  after  each  missile  firing  run,  ISEAS  will 
return  to  the  original  missile  parameters. 

B.  ENCOUNTER  TERMINATION 

1 .  Flight  Termination 

Missile  flight  will  be  terminated  if  any  of  the  following 
conditions  are  encountered: 

a.  Missile  range  exceeds  maximum  range. 

b.  Missile  altitude  exceeds  maximum  altitude. 

c.  Missile  time-of-f light  exceeds  maximum  time  of  flight. 

2 .  Guidance  Termination 

The  ability  to  terminate  missile  guidance  due  to  a  target 
tracker  break-lock  is  not  incorporated . 

C.  MISSILE  POSITION  IN  FLIGHT 

1 •  Global  Coordinate  System 

The  global  cartesian  coordinate  system  is  used  to  define  the 
target  location.  MICE  II  specified  no  reference  point  for  its  global 
coordinate,  so  for  the  ISEAS  missile  module  the  global  coordinate  is 
attached  on  sea  surface  with  the  Y  vector  to  the  north,  the  X  vector  to 
the  east,  and  the  Z  vector  perpendicular  to  X  and  Y  vectors  of  the 
X , Y , Z  coordinates  (See  Figure  3). 

2 .  Missile  Stability  Coordinate  System 

The  missile  stability  axis  system  is  used  to  orient  the  missile 
in  the  global  space,  and  its  origin  is  at  the  missile’s  center  of 
gravity.  The  x  axis  is  along  the  missile  velocity  vector  (with 
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X 


Figure  3.  Coordinate  Systems 
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positive  forward),  the  y  axis  is  along  the  right  wing,  and  the  z  axis 
is  perpendicular  to  the  x  and  y  axis  (with  positive  downward).  See 
Figure  3. 

The  coordinates  X,Y,Z  locate  the  the  missile’s  center  of  gravity 
with  respect  to  the  global  coordinates.  The  angles  gamm  and  etam, 
shown  in  Figure  3,  are  the  elevation  and  azimuth  angles  of  the 
missile’s  velocity  vector  with  respect  to  the  horizontal  plane  of  the 
global  coordinates.  Once  the  missile  orientation  is  defined  in  the 
global  coordinate  system,  transf ormation  to  the  missile  stability  coor¬ 
dinate  system  can  be  performed  using  the  following  direction  cosine 
matrix  in  the  C  language  format.  [Ref.  1:p.  4] 


~T[0][0]  T [0] [ 1 ]  T[ 0][2T" 

= 

T[ 1 ]  [0]  T[1][1]  T[ 1 ] [2] 

-  i3  — 

_T[2][0]  T  [2]  [  1  ]  T[2][2]_ 

__iz_ 

where 


IX’  Iy  Iz  ■ 

ix,iy,ij  = 

T[0]  [0]  = 

T [0]  [1 ]  = 

T [ 0] [2]  = 
T[1][0]  = 

T  [  1  ]  [  1  ]  = 
T [ 1 ] [2]  = 


the  unit  global  coordinate  vectors 
the  unit  missile  coordinate  vectors 

cos ( gamm) *cos( etam) 
cos( gamm) *s in (etam) 
sin( gamm) 
sin( etam) 
cos(etam) 

0 


(4) 
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T [ 2 ] [ 0 ]  =  -sin(gamm)*cos(etam) 

T[2][l]  =  -sin(gamm)*sin(etam) 

T[2][2]  =  cos(gamm) 

For  example, 

ix=  cos ( gomm )#cos( etam)# 1^  +  cos( gamm)*sin( etam)*Iy+  sin(gamm)*l£ 
is  the  unit  vector  of  the  missile  stability  coordinate  in  the  direction 
of  the  missile  velocity  vector. 

D.  LAUNCHER  PARAMETERS 

A  missile  launch  requires  the  following  conditions: 

1.  Line  of  sight  range  (from  the  missile  to  target)  must  be  less  than 
or  equal  to  the  maximum  target  lock-on  range  of  seeker. 

2.  Line  of  sight  range  must  be  less  than  the  maximum  missile  launch 
range . 

3.  Distance  to  the  target  after  total  missile  flight  time  is  less 
than  or  equal  to  the  maximum  missile  range. 

MICE  II  missile  launchers  are  stationary  and  have  no  launcher  bearing 

constraints.  Since  the  ISEAS  module’s  missile  launcher  is  located 

aboard  ship,  and  realistic  shipboard  launch  conditions  are  desired,  the 

following  limitations  were  added  to  the  missile  launcher  requi rements . 

4.  Launcher  bearings  must  not  exceed  135  degrees  port  or  starboard  of 
ship’s  head. 

5.  Launcher  elevation  must  be  less  than  or  equal  to  90  degrees.  This 
will  permit  a  vertical  launch  if  desired. 

If  any  of  the  above  requirements  are  not  satisfied,  missile  launch¬ 
ing  will  not  be  permitted,  and  target  tracking  and  intercept  solution 
calculations  will  continue  until  all  conditions  are  met.  If  launch 
conditions  are  satisfied,  missile  launch  will  be  permitted. 
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The  m_lchr  function  computes  the  launch  angles  (Azi  and  Elev)  based 
on  current  ship  and  calculated  target  intercept  positions.  See  Figure 
4  for  the  m_lchr  block  diagram.  Since  the  launcher  is  located  on  the 
ship,  and  the  ship  is  in  motion,  launch  angles  will  be  relative  to 
ship’s  heading.  The  launch  angles  computation  is  interpreted  from  MICE 
II  but  with  added  launcher  constraints  as  discussed  above.  The  launch 
angle  equations  are  [Ref.  1:p.  11] 

YTINT  -  YMSL 

Azi  =  Arc  tan  [ - ]  (5) 

XTINT  -  XMSL 

ZTINT  -  ZMSL 

Elev  =  Arctan  [ - ]  (6) 

sqrt ( (  XTINT  -  XMSL  )  +  (  YTINT  -  YMSL  )  ) 

where 

XMSL , YMSL , ZMSL  =  the  missile’s  initial  position 

XTINT , YTINT , ZTINT  =  the  target’s  position  at  time  of 

intercept  (provided  by  others). 

If  the  azimuth  or  elevation  limits  are  exceeded,  m_lchr  will  return 
to  the  main  function  in  ISEAS  for  recomputation  of  intercept  solutions 
or  possibly  for  the  user  to  maneuver  the  ship  for  better  launch  angles. 
If  Azi  and  Elev  are  within  the  azimuth  and  elevation  limits,  the  pro¬ 
gram  will  next  check  if  the  target  is  within  range.  If  out  of  range, 
m_lchr  will  return  to  the  main  function  and  the  radar  will  continue 
tracking  the  target.  If  the  target  is  within  range,  m_lchr  will  return 
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Figure  4.  M_Lchr  Flow  Diagram 
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with  the  calculated  azimuth  and  elevation  values.  This  will  be 
converted  to  angles  in  the  global  coordinate  system  and  will  be  the 
initial  gamm  and  etam  of  the  missile. 

E.  MISSILE  FLYOUT  PHASE 

There  are  four  functions  that  m_flyout  calls  to  compute  the 
missile’s  linear  acceleration  and  angular  velocities.  These  functions 
are  the  m_atm(),  m_cna(),  m_thrust(),  and  m_guide().  M_flyout  is 
discussed  after  these  four  functions. 

1 .  M  Atm( ) 

The  m_atm  function  determines  the  atmospheric  densities  and 
temperatures  as  a  function  of  missile  altitude  (ZMSL).  MICE  II  employs 
a  table  look-up  routine  to  determine  the  atmospheric  temperature  and 
density.  The  temperatures  and  air  densities  in  ISEAS  are  approximated 
by  the  following  linear  regression  equations  [Ref.  3:p.  290]  based  on 
the  atmospheric  data  of  Table  1  [Ref  1:p.  47].  See  Figure  5  for  m_atm 
flow  diagram. 


temp  = 

287.25  - 

0 . 00 624 185* ZMSL ,  for 

ZMSL  < 

11500  m 

temp  = 

216.6, 

for 

ZMSL  > 

1 1500  m 

rho  = 

1 .0177572 

-  0 . 000048539*ZMSL 

where 

temp  =  temperature  of  air  at  altitude  ZMSL  [  Kelvin  ] 
rho  =  density  of  air  at  altitude  ZMSL  [  kg/cu.  m  ] 
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TABLE  1 


ATMOSPHERIC  DATA 


ALTITUDE 

(FT) 

-2 000. CO 
0.00 
3000.00 
<1000.00 
6000.00 
5000. 00 
10000 • 00 
12000.00 
14030.00 
160C0. 00 

laccc. oo 
20000. 00 
2Z003. OQ 
24000.00 
26000  .  0  0 
23000. 00 
30003. 00 
J  5  0  0  0  .  0  0 
40000 • 30 
45000.00 
50000.00 
55000. CO 
60300. 00 
65303 .00 
7000  0 . 3  0 
7  5  C 0  0  .  00 
60C00. 00 
65C0Q. 00 
90000. 00 
95000. 00 
103300.00 
1 C50C0 . 30 
1100C3.30 
1  15C03. 30 
120330.00 
125 3 3C. 30 

130003.30 
135300.00 
14CCC0. CO 
145000. 00 

150000.30 
155C3G.03 
» 6CC3C.OO 
165030.00 
1TC3CC.OO 
175CCC.30 
133300.00 
13500 J . 30 
19 J033. 1 0 
195 JOO.OO 
200300.00 
205300.00 


TEhp£*ATU?€ 
(DZ3  K) 


292 

.12 

200 

•  16 

234 

.20 

230 

.24 

276 

.21 

272 

•  32 

260 

.  36 

264 
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Figure  5.  M_Atm  Flow  Diagram 
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2.  M  Cna (  ) 


The  m_cna  function  computes  the  slope  of  the  curve  that  relates 
the  lift  on  the  missile  to  the  angle  of  attack  (the  lift  curve  slope, 


cna)  for  all  missiles 

[Ref.  2 : p .  295].  MICE  II  used  a  table  look-up 

for  cna.  In  ISEAS,  a 

different  method  is  used. 

The  lift  curve  slope  is  dependent  upon  the  missile  altitude 
The  sonic  velocity  at  altitude  is  computed  using 

accel  =  sqrt(  gamma  *  g  *  R  *  temp  )  (9) 


where 

gamma 

=  1.4 

9 

=  gravity  const.  [  kg-m/N-sq.  sec  ] 

R 

=  287  [  N-m/kg-Kelvin  ] 

The  missile  mach  number  is  determined  by 


mach 

=  vmsl  /  accel  (10) 

where 


vmsl 

=  missile  velocity  [  m/sec  ] 

The  lift  curve  slope  coefficient  is  computed  using  the  equation 


cna  * 

4  /  b  (11) 

where 

b  =  sqrt(  mach  *  mach  -  1  )  (12) 

See  Figure  6  for  m_cna  flow  diagram. 
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Figure  6.  M_Cna  Flow  Diagram 
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3 .  M  Guide(  ) 


The  m_guide  function  computes  the  missile  command  angles,  alpha 
(missile  vertical  angle  of  attack)  and  beta  (missile  side  slip  angle). 
See  Figure  7  for  the  m_guide  flow  diagram. 

The  MICE  II  Proportional  Navigation  guidance  equations  are  used 
in  computing  the  angles,  alpha  and  beta,  which  direct  the  missile’s 
flight  path.  The  solutions  for  alpha  and  beta  are  base  upon  the  dif¬ 
ference  in  the  target  position  between  two  consecutive  time  steps.  For 
example,  the  closing  velocity  (closv)  is  determined  by  taking  the 
difference  of  the  range  between  the  missile  and  the  target  from  the 
previous  time  step,  to  the  range  of  the  current  time  step.  Thus, 
closv  =  (  orng  -  rng  )  /  dt  (13) 


where 

closv  *  closing  velocity  [m/sec] 
rng  =  new  missile  to  target  range  [m] 
orng  =  old  missile  to  target  range  [m] 
dt  =  time  increment  [sec] 


The  direction  cosines  of  the  line  of  sight  vector  of  the  new  time  step 
in  the  global  coordinate  system  were  computed  as  follows: 


II 

1 — 1 

i _ i 

> 

XTAR 

-  XMSL 

)/ 

rng 

in 

X 

direction 

V[1]  =  ( 

YTAR 

-  YMSL 

)/ 

rng 

in 

Y 

direction 

V[2]  =  ( 

ZTAR 

-  ZMSL 

)/ 

rng 

in 

Z 

direction 

where 

XT AR , YTAR , ZTAR  =  target  position 
XMSL , YMSL , ZMSL  =  missile  position 
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Figure  7.  M_Guide  Flow  Diagram 
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The  rate  of  change  of  the  line  of  sight  is  computed  using 
losdt[n]  =  (  V[n]  -  oV[n]  )  /  dt 


(15) 


where 

losdt[n]  =  rate  of  change  of  the  line  of  sight 
for  the  nth  direction 

V[n]  =  new  direction  cosine  of  line  of  sight 
oV[n]  =  old  direction  cosine  of  line  of  sight 
n  =  0,1,2 

Using  the  tranf ormation  matrix  presented  in  Section  C,  the  rate  of 
change  of  the  line  of  sight  is  transformed  from  the  global  coordinate 
system  to  the  missile  stability  coordinate  system. 

m[i]  =  m[i]  +  t[i][j]  *  losdt[j]  ,  j  =  0,1,2  (16) 

where 

m[i]  =  rate  of  change  of  the  line  of  sight 
along  the  ith  direction 

t[i][j]  =  transf ormation  matrix  elements 
provided  in  Section  C 

i  =  0,1,2  0  =  x  direction 

1  =  y  direction 

2  =  z  direction 


According  to  MICE  II,  alpha  is  given  by  [Ref.  1:p.  20] 

(  k*closv*m[2]  )  *  W 

alpha  - -  (17) 

cna  *  (S/2)  *  rho  *  vmsl  *  vmsl 


where 


m[2]  =  time  rate  of  change  of  the  direction  cosine 
of  the  LOS  vector  in  the  z  axis 

W  =  missile  weight  [kg] 

S  =  missile  reference  area  [  sq .  m  ] 
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Similarly,  beta  is  determined  by, 


(  k  *  closv  *  m[l]  )  *  W 

beta  - -  (18) 

cna  *  (S/2)  *  rho  *  vmsl  *  vmsl 


where 

m[l]  =  time  rate  of  change  of  the  direction  cosine 
of  the  LOS  vector  in  the  y  axis 

5.  M  Thrust( ) 

The  m_thrust  function  determines  the  missile  thrust  and  missile 
weight  as  a  function  of  time.  NICE  II  used  various  numerical  constants 
on  its  missile  thrust  and  weight  computations  that  were  not  thoroughly 
discussed.  In  ISEAS,  missile  thrust  is  taken  as  a  constant  value  for 
each  type  of  missile  selected,  and  zero  when  the  burn  time  is  exceeded. 
The  missile  weight  is  computed  using 


W  =  oW  -  (  burn  rate  *  dt  )  (19) 

where 

W  =  new  missile  weight  [  kg  ] 

oW  =  old  missile  weight  [  kg  ] 

burn  rate  =  missile  burn  rate  (depends  on  selected 
missile  type)  [  kg/sec  ] 

dt  =  missile  flight  time  increment 

See  Figure  8  for  the  m_thrust  flow  diagram. 

6.  N  Flyout (  ) 

Once  the  missile  is  fired  by  the  user,  the  m_flyout 
function  will  check  the  following  missile  flight  conditions  at 
each  time  step: 
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Figure  8.  M_Thrust  Flow  Diagram 
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1 .  Missile  altitude  will  be  compared  to  the  maximum  missile  altitude 
of  20 , 000  meters . 

2.  Time  of  flight  will  be  compared  to  the  maximum  missile  flight 
time . 

3.  Range  between  the  missile  and  the  target  will  be  computed  and 
compared  to  the  maximum  separation  range  allowed. 

If  any  of  the  above  flight  checks  are  exceeded,  the  rn_flyout  function 

will  stop  and  return  to  the  main  program. 

Another  check  is  the  minimun  range  between  the  missile  and  the 
target.  If  the  range  is  6.1  meters  or  less,  this  will  be  considered  a 
hit.  See  Figure  9  for  the  m_flyout  flow  diagram. 

Once  alpha  and  beta  are  computed  by  the  m_guide  function,  lift 
conditions  are  generated  and  the  linear  acceleration  equation  for  the 
vdotm  is  [Ref.  1:p.  6] 


[T  -  D  -  W  *  sin(gamm)  -  La  *  alpha  -  Lb  *  beta  ] 

vdotm  =  - 

W  /  g 


(20) 


where 

T  =  missile  thrust  [N] 

D  =  missile  drag  at  zero  lift  [N] 

La  =  cna  *  q  *  S  *  alpha 

missile  lift  due  to  angle  alpha 

Lb  =  cna  *  q  *  S  *  beta 

missile  lift  due  to  angle  beta 

q  =  dynamic  pressure  =  . 5*rho*vmsl*vmsl  [N/sq.  m] 


See  Figure  10  for  the  forces  on  the  missile. 
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Figure  9.  M_Flyout  Flow  Diagram 
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Calculation  of  the  missile  elevation  angle  turn  rate  (gamdm) 


due  to  the  angle  of  attack,  alpha,  gravity  and  the  lift  component  due 
to  the  net  thrust  is  determine  by  using  [Ref.  1:p.  6] 


[  ( T  -  D )  *  alpha  +  La  -  W  *  g  *  cos(gamm)  ] 

gamdm  =  -  (21) 

W  *  vmsl 


Similarly  the  missile  azimuth  angle  turn  rate  (etadm)  due  to 
the  side  slip  angle,  beta,  is  given  by  [Ref.  1:p.  6] 


etadm 


[  -(T  -  D )  *  beta  -  Lb  *  beta  ] 


W  *  vmsl  *  cos(gamm) 


(22) 


Figure  10.  Forces  On  Missile 
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From  the  computed  linear  acceleration  and  angular  velocities, 
an  update  of  missile  position,  orientation  and  velocity  is  determined. 
In  MICE  II  [Ref.  1 : p p .  4-10],  a  slope  averaging  technique  is  used  to 
numerically  update  the  missile  velocity  and  location. 

The  missile  is  at  initial  position  at  time  of  fire.  Thus, 


XMSL  =  XC 

vmsl  =  VC 

YMSL  =  YC 

ogamm  =  gamm 

ZMSL  =  ZC 

oetam  =  etam 

where 


XC, YC.ZC 
VC 
ogamm 
oetam 


values  used 
values  used 
old  missile 
old  missile 


in  position  corrections 
in  velocity  corrections 
elevation  angle 
azimuth  angle 


At  the  first  time  step  (dt),  the  missile  is  in  flight  and  com¬ 
putation  of  the  linear  acceleration  and  angular  velocities  provide 
updates  to  missile  position,  orientation  and  velocity.  [Ref.  1:p.  9] 


vmsl  =  VC  +  dt 

*  vdotm 

(23) 

gamm  =  ogamm  + 

dt  * 

gamdm 

(24) 

etam  =  oetam  + 

dt  * 

etadm 

(25) 

XMSL  =  XC  +  .5 

*  dt 

*  vmsl  *  cos(gamm) 

*  cos(etam) 

(26) 

YMSL  =  YC  +  .5 

*  dt 

*  vmsl  *  cos(gamm) 

*  sin(etam) 

(27) 

ZMSL  =  ZC  +  .5 

*  dt 

*  vmsl  *  sin(gamm) 

(28) 
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The  following  are  saved  after  the  computation  of  equations  (23) 
to  (28):  XC,  VC,  ZC ,  gamdm,  etadm,  vmsl,  etam,  gamm,  vdotm,  then  are 
used  as  the  old  values,  oXC,  oYC,  oZC,  ogamdm,  oetadm,  ovmsl,  oetam, 
ogamm,  ovdotm,  for  the  next  time  step. 

As  the  missile  continue  on  its  flight  path,  at  time  greater 
than  zero,  gamdm  and  etadm  values  are  generated  and  corrections  are 
computed  by  the  slope  average  technique.  [Ref.  1:p.  9] 


GAMC  =  ogamm  +  . 5  *  dt  *  (ogamdm  +  gamdm) 

ETAC  =  oetam  +  .5  *  dt  *  (oetadm  +  etadm) 

VC  =  oVC  +  . 5  *  dt  *  (ovdotm  +  vdotm) 

XC  =  oXC  +  .5  *  dt  *  (oVC  +  VC)  *  cos(GAMC)  *  cos(ETAC) 

YC  =  oYC  +  .5  *  dt  *  (oVC  +  VC)  *  cos(GAMC)  *  sin ( ETAC ) 

ZC  =  oZC  +  .5  *  dt  *  (oVC  +  VC)  *  sin(GAMC) 


(29) 

(30) 

(31) 

(32) 

(33) 

(34) 


where 


GAMC  =  correction  value  for  elevation  angle 
ETAC  =  correction  value  for  azimuth  angle 


The  missile  position, 


orientation  and  velocity  is  updated  by, 


[Ref.  1 : p .  10] 


vmsl 

gamm 

etam 

XMSL 

YMSL 

ZMSL 


oVC  +  d t  *  vdotn 

ogamm  +  dt  *  gamdm 

oetam  +  dt  *  etadm 

XC  +  .5  *  dt  *  (VC  +  vmsl) 

YC  +  . 5  *  dt  *  ( VC  +  vmsl) 

ZC  +  .5  *  dt  *  (VC  +  vmsl) 


*  cos(gamm) 

*  cos(gamm) 

*  sin(gamm) 


(35) 

(36) 

(37) 

*  cos(etam)  (38) 

*  sin(etam)  (39) 

(40) 
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The  computed  XC,  YC,  ZC,  ETAC,  GAMC ,  vdotm,  vmsl,  etadm,  and  gamdm 


are  saved  and  becomes  the  old  values  for  the  next  time  step. 

The  output  from  this  function  will  consist  of  missile  position 
( XMSL , YMSL , ZMSL ) ,  elevation  angle  (gamm),  azimuth  angle  (etam),  velo¬ 
city  (vmsl),  weight  (W),  and  missile  thrust  (T).  All  of  which  the  main 
function  can  utilize  for  graphical  presentations. 

F.  KILL  PROBABILITY 

This  portion  will  conduct  the  endgame.  It  will  be  developed  in  a 
future  study. 
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IV.  CONCLUSIONS  AND  RECOMMENDATIONS 


The  C  programing  language  was  used  for  this  model  because  of  its 
advertised  user  friendly  dialogue  techniques  and  its  simple  and  effi¬ 
cient  operations.  Although  the  simplicity  and  generality  favor  the  C 
programing  language,  since  there  is  no  standard  definitions  for  the  C 
language,  compiling  using  Lattice  C  and  Microsoft  C  compilers  causes 
errors  found  from  one  unit  but  not  from  the  other.  It  is  therefore 
recommended  to  work  with  one  selected  compiler  until  the  language  is 
standardized . 

MICE  II  simulation  logics  was  easily  followed  and  understood.  The 
translation  to  the  C  language  for  the  thrust  and  weight  calculations 
presented  some  problems  because  of  the  numerous  undocumented  constant 
variables  used  in  the  equations.  The  use  of  the  table  look-up  method 
also  created  difficulties  in  translating,  therefore,  alternate  equa¬ 
tions  were  used  that  provided  basically  the  same  results.  All  of  the 
other  major  equations  provided  no  additional  difficulties  in  C  language 
translation . 

The  use  of  the  missile  model  in  decision  making  can  assist  the  user 
by  providing  trade-offs  and  relationships  and  by  evaluatiing  and  com¬ 
paring  alternatives.  But  further  work  remains  to  be  done.  The  missile 
model  is  currently  working  with  only  one  guidance  mode.  Hopefully,  ad¬ 
ditional  guidance  modes  (command  to  line  of  sight  guidance,  beam  rider, 
and  track-via-missile )  can  be  implemented. 
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Finally,  during  the  early  stages  of  this  project,  it  is  recommended 


that,  (1)  the  student  develope  a  working  understanding 
of  the  equations  of  motion  and  guidance  law,  and  (2) 
complete  a  formal  C  language  programing  course. 


on  the  mechanics 
that  the  student 


35 


LIST  OF  REFERENCES 


1.  Chan,  P.T.  and  Huffman,  R.A.,  Surface  To  Air  Missile 
Model  -  Mice  II ,  Vol  I,  Analyst  Manual,  Vought  Corporation, 
Dallas,  Texas,  April  1980. 

2.  Lindsey,  Gerald  H.  and  Redmon ,  Dan  R. ,  Tactical  Missile  Design , 
notes  presented  for  AE4703  Missile  Stability  and  Design,  Naval 
Postgraduate  School,  Monterey,  California,  July  1986. 

3.  Miller,  I.  and  Freund,  J,E.,  Probability  And  Statistics  For 
Engineers ,  3rd.  Edition,  Prentis-Hall ,  Inc.,  Inglewood  Cliffs, 
New  Jersey,  1985. 


36 


APPENDIX  A 


The  following 

alphabetical  order 

accel  [m/sq.sec] 
alpha  [rad] 
altmax  [meters] 
azi  [rad] 

beta  [rad] 

cd 

cl 

closv  [m/sec] 
cna 

dt  [sec] 
drag  [N] 

elev  [rad] 
etam  [rad] 
etac 

etadm  [rad/sec] 

g[kg-m/N-sq.sec] 

game 

gamm  [rad] 
gamdm  [rad/sec] 

ko 

lift  [N] 
losdt  [rad/sec] 

m[]  [m] 

mach 

mbrl  [kg/sec] 
mbrm  [kg/sec] 
mbrs  [kg/sec] 
m_l  r 
m_mr 
mselect 
m  sr 


PROGRAM  VARIABLES 

is  a  listing  of  computer  program  variables  in 


missile  acceleration 

missile  attitude  angle  wrt  stability  axis 
maximum  altitude 

relative  azimuth  angle  of  missile  launcher 

missile  side  slip  angle  wrt  stability  axis 

missile  drag  coeff 
missile  lift  coeff 
missile-target  closing  velocity 
missile  lift  curve  slope  coeff 

time  increment 
missile  drag 

relative  elevation  angle  of  missile  launcher 
azimuth  angle  of  missile  velocity  vector 
correction  value  for  azimuth  angle 
rate  of  change  of  missile  azimuth  angle 

gravity 

correction  value  for  elevation  angle 
elevation  angle  of  missile  velocity  vector 
rate  of  change  of  missile  elevation  angle 

base  value  of  Proportional  Navigation  constant 


missile  lift 

rate  of  change  of  LOS 


missile  position  in  stability  coord 

missile  mach  number 

long  range  missile  burn  rate 

med  range  missile  burn  rate 

short  range  missile  burn  rate 

long  range  missile  (3) 

med  range  missile  (2) 

missile  select  by  user 

short  range  missile 
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oetadm 

oetam 

ogamdm 

ogamm 

orng  [meters] 


old  rate  of  change  azimuth  angle 
old  azimuth  angle 

old  rate  of  change  elevation  angle 
old  elevation  angle 

old  range  between  missile  and  target 


ov[]  [rad] 

old 

direction  < 

cosine 

of 

LOS 

ovc 

old 

correct ion 

value 

for 

miss 

oxc 

old 

correction 

value 

for 

xmsl 

oyc 

old 

correction 

value 

for 

ymsl 

ozc 

old 

correction 

value 

for 

zmsl 

rho  [kg/cu.m] 
rng  [meters] 
rngl  [m] 
rng2  [m] 
rng3  [m] 
rtint  [m] 

s  [sq .  meters] 


dynamic  pressure 

atmospheric  density 
range  between  missile  and  target 
max  range  of  short  range  missile 
max  range  of  med  range  missile 
max  range  of  long  range  missile 
target  intercept  range 

surface  area  of  missile 


t  [secs] 
tl  [secs] 
t2  [secs] 
t3  [secs] 
temp  [kelvin] 


missile  flight  time 

max  flight  time  for  short  range  miss 
max  flight  time  for  med  range  missil 
max  flight  time  for  long  range  missi 
atmospheric  temperature 


thru 

st  [N] 

missile 

th 

rust 

thru 

stl 

[N] 

thrust 

of 

short  i 

range  missile 

thru 

st2 

[N] 

thrust 

of 

med  rai 

nge 

missile 

thru 

st3 

[N] 

thrust 

of 

long  range 

missile 

v[] 

[ra 

d] 

directi 

on 

cosine 

of 

LOS 

VC 

correct 

ion 

i  value 

for 

vmsl 

vmsl 

[m/ 

s] 

missile 

velocity 

ile 

e 

le 


w  [kg] 
wlr  [kg] 
wmr  [kg] 
wsr  [kg] 


missile  weight 
weight  of  long  range  miss 
weight  of  med  range  missi 
weight  of  short  range  mis 


ile 

le 

sile 


xc 

xmsl  [m] 
xtar  [m] 
xtint  [m] 


correction  value  for  xmsl 

missile  position  on  x  axis 

target  position  on  x  axis 

target  intercept  position  on  the  x  axis 


yc 

ymsl  [m] 
ytar  [m] 
ytint  [m] 


correction  value  for  ymsl 

missile  position  on  y  axis 

target  position  on  y  axis 

target  intercept  position  on  the  y  axis 
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zc 

zmsl  [m] 
ztar  [m] 
ztint  [m] 


correction  value  for  zmsl 

missile  position  on  z  axis  (missile  altitude) 

target  position  on  z  axis 

target  intercept  position  on  the  z  axis 
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APPENDIX  B 


PROGRAM  LISTINGS 

The  following  is  a  "C"  language  programming  listing  of  the 
ISEAS  missile  simulation  program. 
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m_lchr( ) 

/*  This  function  will  compute  the  launcher’s  azimuth  and  elevation 
angles  with  respect  to  the  ship’s  heading.  Input  to  this 
function  will  the  launcher  (the  initial  missile  position) 
position  and  the  calculated  target  intercept  position.  The 
output  will  be  the  launcher’s  azimuth  and  elevation  angles 
and  must  be  converted  to  the  global  coordinates  and  becomes 
the  initial  gamm  and  etam  for  m_flyout.  */ 


{ 

/*  COORDINATES  ARE  RELATIVE  TO  SHIP'S  HEAD  */ 

float  xtin t , ytint , ztint ;  /*tgt  intercept  position*/ 
float  xtar , ytar , ztar ;  /*tgt  rel  position*/ 

float  xmsl , ymsl , zmsl ;  /^launcher  or  initial  missile  position*/ 
float  elmts=90.;  /*max  lehr  elevation*/ 

float  almts=135.;  /*lchr  port/stbd  bearing  limits*/ 
float  rng;  /*  range  to  target  */ 

float  rtint  /*  intercept  range  */ 

float  rngl =20000 , rng2=48000 , rng3= 1 2 1 000 ;  /*  missile  ranges  */ 
float  elev,  azi; 
float  x , y , z , xy ; 

/****  COMPUTE  LAUNCHER  BEARING  ANGLE  *********/ 
azi  =  atan ( ( ytint-ymsl ) / ( xtin t-xmsl ) ) ; 
azi  =  ( 1 80*azi ) / 3 . 1 4 1 59 ; 
if  (abs(azi)  >  almts) 

{ 

printf( " \n  BEARING  ANGLE  LIMITS  EXCEEDED”); 
retu  rn ; 

) 


/****  COMPUTE  LAUNCHER  ELEVATION  ANGLE  ******/ 
x  =  xtint-xsml; 
y  =  ytint-ysml; 
z  =  ztint-zsml 
xy  =  sqrt ( x*x  +  y*y); 
elev  =  atan(z/xy); 
elev  =  ( 1 80*elev ) /3 . 1 41 59  ; 
if  (elev  >  elmts) 

{ 

printf( " \n  ELEVATION  LIMITS  EXCEEDED”); 
return ; 

) 


/****  CHECK  IF  WITHIN  MISSILE  RANGE  ******/ 
rtint  =  sqrt(x*x  +  y*y  +  z*z); 
if  (  mselect  ==  1  ) 

{ 
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if  (  rtint  >  rngl  ) 

{ 

printf( "\nTARGET  NOT  WITHIN  SEEKER  RANGE"); 
printf ( "\nCONTINUE  TRACKING"); 
return ; 

) 

) 

if  (  mselect  =  =  2  ) 

{ 

if  (  rtint  >  rng2  ) 

{ 

printf( " \nTARGET  NOT  WITHIN  SEEKER  RANGE"); 
printf ( " \nCONTINUE  TRACKING"  )  ; 
return ; 

) 

> 

if  (  mselect  ==  3  ) 

( 

if  (  rtint  >  rng3  ) 

{ 

printf( "\nTARGET  NOT  WITHIN  SEEKER  RANGE"); 
printf ( "\nCONTINUE  TRACKING"); 
return ; 

> 

) 


/***»  PRINT  LAUNCHER  AZIMUTH  AND  ELEVATION  »***/ 
if  (  azi  <  0  ) 

printf  (  " \n  LAUNCHER  BEARING  =  f6.2f  DEG  TO  PORT", azi); 
else 

printf("\n  LAUNCHER  BEARING  =  %.2f  DEG  TO  STBD" , azi ) ; 
printf ( "\n  LAUNCHER  ELEV  =  $.2f  DEG" , elev )  ; 
printf ("\n  TARGET  WITHIN  RANGE"); 
pr intf ( " \n  LAUNCHER  READY"); 


) 
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#include<stdio.h> 
#include<math . h  > 


/*  This  program  was  created  to  simulate  the  main  function  and 
provide  target  data  to  evaluate  variable  inputs  and  outputs 
of  the  flyout  function.  Other  functions  are  also  called  upon 
to  provide  the  necessary  values  for  the  computations  of  the 
missile’s  equations  of  motion  and  eventually,  the  missile’s 
position  in  global  coordinates.  */ 


*  GLOBAL  VARIABLES  */ 

float  xtar= 1 0000 , ytar= 1 0000 , ztar= 1 0000 ;/ ^current  target  coord*/ 

float  xmsl , ymsl , zmsl ; /^current  missile  coord  */ 

float  vmsl=671 ;  /*  missile  velocity  =  2200  ft/sec  */ 


float  gamm, ogamm; 
float  gamdm, ogamdm ; 
float  etadm , oetadm ; 
float  etam,oetam; 
float  w; 
float  s; 
float  s 1 =  2 . 3 3 ; 
float  s2=4.92; 
float  s3=8 . 78 ; 
float  thrust; 


/*  elev  angle  of  vmsl  */ 


/' 


'/ 


float  drag , cd  =  0 . 2 ; 

/ 


/*  azi  angle  of  vmsl  */ 
msl  weight  */ 

/*  area  of  selected  msl  */ 

/*  ref  area  of  short  rng  msl 
/*  ref  area  of  med  rng  msl  */ 

*  ref  area  of  long  rng  msl  */ 

*  msl  thrust  [  N  ]  */ 

/*  msl  drag, drag  coeff  */ 

*  msl  lift, lift  coeff  */ 

/*  time  increment  delta-t  */ 

/*  msl  maximum  altitude  in  meters  */ 
msl  flight  time  */ 


float  lif t , cl ; 
float  dt; 

float  altmax=21 000 ; 
float  t;  /* 

float  rng; 

float  orng=20000;  /*  initial  value  from  main  function  */ 

float  rho , temp , g= 1 ; 
float  alpha, beta; 
float  cna; 

float  ko;  /*  input  base  value  of  proportional  constant*/ 

int  mselect; 

float  v  [3] , ov [3] ; 

float  cpa; 

int  shoot; 

int  ans , ans 1 , ans2 , ans3 ; 

float  vc,xc,yc,zc; 

float  ovc , oxc , oyc , ozc  ; 

float  vdotm , ovdotm; 

float  game , etac , ogamc , oetac  ; 

float  wsr=23 1 . 3 , wmr=581 ,wlr=1066; 

float  mbrs= 1 . 36 , mbrm=0 . 338 , mbrl=0 . 338 ; 

float  tl =30 , t2=71 , t3=1 81 ; 

float  thrustl =200000, thrust2=400000,thrust3=400000; 


43 


main(  ) 

{ 

printf("\n  Change  target  position?  1=yes/0  =  no  ”  ) ; 
scant  (  "id  ” ,  &ans); 
if  (  ans  ==  1  ) 

{ 

pr intf ( " \nEnter  target  position  in  meters”); 

printf(”\n  xtar  =  ” ) ; 

scanf(”$f”,  &xtar); 

printf(”\n  ytar  =  ”  ) ; 

scanf(”^f" ,  &ytar); 

printf(”\n  ztar  =  "); 

scant ( "if" ,  &ztar  ) ; 


} 

printf ( ” \n 
printf( " \n 
printf(”\n 
printf(” \n 


Initial 

ta 

rget 

position 

in 

xtar  = 

if 

ytar 

=  if  zta 

r  = 

Target 

hea 

ding 

is  fixed 

at 

at  10  meters  decrement 

per 

meters" ) ; 

if  ”  ,  xtar  ,  ytar  ,  ztar  ) ; 
270  degs” ); 
time  step”); 


printf(”\n\n  Enter 
printf(”\n  xmsl  = 
scant ( "if" ,  &xmsl ) 
pr intf ( ” \n  ymsl  = 
scanf("^f",  &ymsl) 
printf(”\n  zmsl  = 
scanf(”$f”,  &zmsl) 


missile  initial 
"); 

"); 


position  in  meters”); 


printf("\n  Enter  missile  attitude  angle  in  radians  =  ”); 

scanf(,,^f")  &gamm); 

printf(”\n  Enter  missile  azimuth  angle  in  radians  =  " ) ; 

scant ( "if" ,  &etam) ; 


printf(”\n  Desired  time  increment  in  secs  =  ” ) ; 

scant (  "if  "  ,  &dt ) ; 

printf(”\n  Enter  proportional  constant  base  value  =  ”  ) ; 

scant ( "if" ,  &ko ) ; 

printf (”\n  Select  type  of  missile  ” ) ; 

printf (” \nShort  Range=  1  Med  Range=  2  Long  Range=  3  "); 
printf(”\n  missile  select  =  ” ) ; 

scanf(”$d”,  &mselect); 


printf(”\n  Change  missile  ref  area?  1=yes/0=no  "); 
scant ( "if" ,  &ans1 ) ; 
if  (  ansi  ==  1  ) 

{ 
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if  ( mselect  ==  1 ) 

{ 

printf("\n  Enter:  si  =  "); 

scant (  "iof  "  ,  &s1  ) ; 

} 

if  (mselect  ==  2) 

( 

printf(M\n  Enter:  s2  =  " )  ; 

scant ( "$f" ,  &  s  2 ) ; 

} 

if  (mselect  ==  3) 

{ 

printf("\n  Enter:  s3  =  "); 

scant  (  "gf  "  ,  &s3); 

} 

} 

printf("\n  Change  missile  thrust?  1=yes/0  =  no  "  ) 
scanf("$d" ,  &ans2); 
if  (  ans2  ==  1  ); 

{ 

if  (  mselect  ==  1  ) 

{ 

printf("\n  Enter:  thrustl  =  " ) ; 

scant ( "$f " ,  &thrustl); 

} 

if  (  mselect  ==  2  ) 

( 

printf(M\n  Enter:  thrust2  =  " ) ; 

scant ( "$f " ,  &thrust2 ) ; 

} 

if  (  mselect  ==  3  ) 

( 

printf(,,\n  Enter:  thrust3  =  "); 

scant ( H$f" ,  &thrust3); 

} 

} 

printf(M\n  Change  missile  burn  rate?  1=yes/0=no  " 
scant ( "$d" ,  &ans3); 
if  (  ans3  ==  1  ) 

{ 

if  (  mselect  ==  1  ) 

{ 

printf(n\n  Enter:  burn  rate  1  =  " ) ; 

scant ( "#f" ,  &mbrs); 

} 
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if  (  mselect  ==  2  ) 

{ 

printf("\n  Enter:  burn  rate  2  =  M) 
scant (  11 ,  &mbrm); 

} 

if  (  mselect  ==  3  ) 

{ 

printf(M\n  Enter:  burn  rate  3  =  " ) 

scanf( M^f" ,  &mbrl); 

} 

> 


printf ("\n  When  ready  to  shoot  "); 
spacebar ( ) ; 
t  =  0  ; 

while  (  t  <  31  ) 

{ 

printf ( M\n\n  t  =  $f",t); 

m_f lyout  (  )  ; 

xtar  =  xtar  -  10; 
printf("\n  xtar  =  $f",xtar); 
ytar  =  ytar; 

printf ("  ytar  =  ^f" ,ytar); 
ztar  =  ztar; 

printf (  "  ztar  =  $f",ztar); 
t  =  t  +  dt; 

} 

printf("\n  End  of  flight  " ) ; 

} 


spacebar( ) 

{ 

printf ( "\n$s" , "  press  spacebar"  ) ; 
while  (  getch(  )  !  =  '  ’  ) ; 
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m_flyout( ) 

/*  This  function  computes  the  missile’s  linear  accelerations 
and  angular  velocities.  Inputs  of  missile  and  target 
positions,  and  missile  general  character istics  are  required. 
Computations  of  missile  attitude  and  azimuth,  velocity, 
weight,  thrust,  and  position  are  the  outputs  and  are  updated 
at  each  time  increment  */ 

{ 

float  a , b ; 
float  x , y , z ; 

x  =  xtar  -  xmsl; 
y  =  ytar  -  ymsl; 
z  =  ztar  -  zmsl; 

/**************  range  between  missile  and  target  **************/ 

rng  =  sqrt(  x*x  +  y*y  +  z*z); 
printf("\n  range  =  %f  meters" , rng); 

/**********  altitude  check  of  missile  flight  ******************/ 

if  (  zmsl  >  altmax  ) 

{ 

printf(M\n  max  altitude  exceeded  " ) ; 
return ; 

} 


/ 


*********  QpA  check  0f  missile  flight 


/ 


if  (rng  <=  6.1  )  /*cpa  .0061  km  or  less  */ 

return;  /*  warhead  exploded-  calc  miss  dist  */ 

if  (rng  >  orng ) 

{ 

printf(M\n  missile  pass  CPA"); 
cpa  =  orng; 

printf("\n  CPA  =  #f  m  ",cpa); 
t  =  31  ; 
return ; 

} 


y***-*-**-**.* 

/********* 


call  function  to  calc  temp  and  density  wrt 
altitude  -  return  value  of  temp  and  density 


***.*.*.*.*.*  j 
********* j 


m_atm( ) ; 
printf(" \n 
printf(" \n 
printf(" \n 


At  an  altitude  of  $f",zmsl); 
temperature  =  %fn , temp ) ; 
density  =  ^f", rho); 
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/*********  call  m_cna  function  to  calc  the  lift  curve  slope  **/ 
/*********  coefficient  (  cna  )  with  respect  to  altitude  and  **/ 
/*********  temperature  **/ 

m_cna( ) ; 

printf("\n  cna  =  $fM,cna); 

/********  call  function  to  calc  msl  propulsion  forces  ********/ 
/********  return  with  thrust  and  msl  weight  ** ***************** */ 

m_thrust ( ) ; 


/*******  call  function  to  calc  msl  guidance  -  return  ***********/ 
/*******  alpha  and  beta  -  msl  command  angles  *******************/ 


m_guide( ) ; 

/**********  calc  drag  at  zero  nft  *****************************/ 

printf(M\n  cd  =  $f",cd); 
printf("\n  rho  =  $f",rho); 
printf(”\n  s  =  ^f",s); 
printf(M\n  vmsl  =  ^f", vmsl); 
drag  =  cd* ( s/2 ) *rho*vmsl * vmsl ; 
printf("\n  drag  =  ^f 11 ,  drag) ; 


/**********  At  time  zero,  missile  is  fire 
/**********  misssile  position  is  launcher 
if  (  t  ==  0  ) 

{ 


and  initial**********/ 
position  ************/ 


pr  i 

ntf ( " \n 

At  time  of 

pri 

ntf ( " \n 

xmsl  = 

pri 

ntf ( "\n 

etam= 

VC 

=  vmsl; 

xc 

=  xmsl; 

yc 

=  ymsl; 

/  ** 

zc 

=  zmsl; 

!  ** 

oetam  =  etam; 
ogamm  =  gamm; 
return ; 


fire"  ) ; 
ymsl  =  °fof 
gamm=  °fof 


values  a 
for  the 


zmsl=  °t> f  "  ,  xmsl ,  ymsl ,  zms 
" , etam, gamm) ; 


ssigned  to  old  values  **/ 
next  time  step  **/ 


l); 


} 


/************* 


calc  of  msl 


linear  acceleration  *****************/ 


a  =  cna* . 5*s*rho*vmsl*vmsl*alpha ; 
b  =  cna* . 5*s*rho*vmsl*vmsl*beta ; 
printf("\n  a  =  °fof  b=  %f  "  ,  a  ,  b  ) ; 

vdotm=((thrust-drag-(w*g*sin( gamm ))-(a*alpha+b*beta))/w); 
printf("\n  vdotm  =  <j6f”  ,vdotm); 
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/  *****.*******  calc  of  msl  elavation  angle  turn  rate  **-**********/ 

gamdm =( (thrust-drag ) *alpha+a*w*g-( w^g*cos( gamm ) ))/(w*vmsl); 
printf("\n  gamdm  =  ^f",gamdm); 

/  ************  ca^c  0f  ms]_  azimuth  angle  turn  rate  *************-*/ 

etadm=( -(thrust-drag)*beta-b) / ( w*vmsl * cos ( gamm ) ) ; 
printf("\n  etadm  =  $f", etadm); 


/****  At  this  time,  the  missile  is  in  flight  and  if  the  time  **/ 
/****  step  equals  the  time  of  increment,  then  first  iteration**/ 
/****  of  missile  position,  orientation,  and  velocity  is  ***/ 

/****  computed.  *** / 

if  (  t  ==  dt  ) 

{ 

/*  update  of  msl  vel  */ 
vmsl  =  vc  +  0 . 5*dt * vdotm ; 


/*  new  msl  elevation  angle  */ 
gamm  =  ogamm  +  0 . 5*dt*gamdm ; 

/*  new  msl  azimuth  angle  */ 
etam  =  oetam  +  0 . 5*dt *etadm ; 

/*  update  msl  position  base  on  flight  correction  */ 
xmsl  =  xc  +  0 . 5*vmsl*dt*cos ( gamm ) *cos( etam ) ; 
ymsl  =  yc  +  0 . 5* vmsl* dt*cos ( gamm ) *sin ( etam ) ; 
zmsl  =  zc  +  0 . 5*vmsl*dt*sin( gamm) ; 

/*  change  new  value  to  old  value  for  next  iteration  */ 


ovdotm 

= 

vdotm 

ogamdm 

= 

gamdm 

oetadm 

= 

etadm 

ovc  = 

vmsl ; 

oxc  = 

xc 

; 

oyc  = 

yc 

; 

ozc  = 

zc 

» 

oetam 

= 

etam ; 

ogamm 

= 

gamm ; 

} 
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/****  since  the  time  step  is  greater  than  the  time  increment  **/ 
/****  and  greater  than  zero,  the  missile  is  well  into  its  **/ 

/****  flight  path  and  the  missile  position,  orientation,  and  **/ 
/*•***  velocity  are  computed,  and  are  corrected  by  the  slope  **/ 
/**-*-*  average  technique  **/ 

else 

{ 

game  =  ogamm  +  . 5  *  dt  *  (ogamdm  +  gamdm); 
etac  =  oetam  +  . 5  *  dt  *  (oetadm  +  etadm); 


VC  = 

ovc  +  . 

5  *  dt  *  (ovdotrr 

l  +  vdotm ) ; 

XC  = 

oxc  +  . 

5  *  dt  *  (ovc  + 

vc  )  * 

cos(gamc)  * 

cos( etac ) ; 

yc  = 

oyc  +  . 

5  *  dt  *  (ovc  + 

vc )  * 

cos(gamc)  * 

sin ( etac ) ; 

zc  = 

ozc  +  . 

5  *  dt  *  (ovc  + 

vc )  * 

sin ( game  ) ; 

vmsl 

=  vc  + 

dt  *  vdotm; 

gamm 

=  ogamm  +  dt  *  gamdm; 

etam 

=  oetam  +  dt  *  etadm; 

xmsl 

=  xc  + 

.5  *  dt  ‘  (vc  + 

vmsl ) 

*  cos(gamm) 

*  cos(etam); 

ymsl 

ii 

< 

o 

+ 

.5  *  dt  *  (vc  + 

vms  1 ) 

*  cos(gamm) 

*  sin(etam); 

zmsl 

=  ZC  + 

. 5  *  dt  *  (vc  + 

vmsl ) 

*  sin( etam) : 

ovdotm  =  vdotm; 
ogamdm  =  gamdm; 
oetadm  =  etadm; 
oetam  =  etac; 
ogamm  =  game; 
ovc  =  vc; 
oxc  =  xc; 
oyc  =  yc; 
ozc  =  zc; 


} 


/*********  Print  missile  position,  orientation,  and  velocity  **/ 
I  *********  as  Weii  as  missile  thrust  and  weight  **/ 

printf ( " \nxmsl=  %f  m  ymsl=  %f  m  zmsl=  %f  m  ", xmsl , ymsl , zmsl ) 

m/sec" , vmsl ) ; 
etam  =  %f  " , gamm , etam) ; 
thrust  =  N" , w, thrust ) ; 


printf("\n  vmsl  = 
printf("\n  gamm  =  %f 


printf(M\n  w  =  %f  kg 
return ; 


} 
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m_atm( ) 

/*  This  function  will  calculate  the  atmospheric  density  and 
temperature  with  respect  to  altitude.  Equations  used  were 
derived  by  the  linear  regression  method.  Maximum  missile 
altitude  is  21000  meters,  and  the  flyout  subroutine  will 
check  if  the  missile  exceeds  maximum  altitude.  Therefore  as 
long  as  the  missile  is  at  21000  meters  or  below,  atmospheric 
density  and  temperature  will  be  computed.  */ 


{ 


y  *****  * 


rho 


density  computation  [  kg  /  cu .  m  ]  **************** * / 

=  (  1.0177572  -  (0.000048539  *  zmsl)  ); 


/******  temperature  computation 


[  kelvin  ] 


*****************  j 


if  (  zmsl 
temp 
else 
temp 


=  11 

500  ) 

!  *  * 

at 

216  . 

66; 

j  *  * 

is 

287. 

25  -  (0. 

. 006241 

this  altitude,  temperature 
constant  at  216.66  k 

865  *  zmsl ) ; 


*  *  j 

*  *  j 


/**  Return  to 
density 


flyout  subroutine  with 


values  of  temperature  and 

**  j 


> 
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m_cna(  ) 

/*  This  function  calculates  the  lift  curve  slope 

coefficeint  with  respect  to  altitude  and  temperature.  CNA 
is  used  for  the  computation  of  alpha,  beta,  and  the  equations 
of  motions.  */ 

{ 

float  accel; 
float  mach; 
float  b; 

/*■***  calc  missile  acceleration  -***-*/ 
accel  =  20.05  *  sqrt(temp); 


/****  calc  missile  mach  number  ****/ 
mach  =  vmsl  /  accel; 

b  =  sqrt( (mach*mach )  -  1); 


/****  calc  lift  curve  slope  coeff  ****/ 
cna  =  4  /  b; 
return ; 


} 
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m_thrust( ) 

/*  The  m_thrust  function  computes  the  missile  weight  during 
flight.  Missile  thrust  is  a  constant  value  and  this  value 
will  depend  on  the  type  of  missile  selected.  At  burn  out, 
missile  thrust  will  be  zero  and  flight  time  is  exceeded. 

Output  from  this  function  will  be  missile  weight  and  missile 
thrust.  *7 

{ 


/***  if  user  selects  a  short  range  missile  ***/ 

if  (  mselect  ==  1  ) 

{ 

s  =  s  1  ; 

if  (  t  ==  0  )  /**  at  missile  time  of  fire  **/ 

{ 

w  =  wsr; 

thrust  =  thrustl ; 
return ; 

> 

else 

{ 

if  (  t  <  tl  )  /**  missile  flight  time  **/ 

{ 

w  =  w  -  (mbrs  *  d t ) ; 
thrust  =  thrustl ; 
return ; 

> 

else 

{ 

printf("\n  missile  burn  time  compeleted" ) ; 
w  =  w; 
thrust  =  0; 
return ; 

) 

> 

) 

/***  user  selects  a  medium  range  missile  ***/ 

if  (  mselect  ==  2  ) 

{ 

s  =  s2 ; 

if  (  t  ==  0  )  /**  missile  time  of  fire  **/ 

{ 

w  =  wmr ; 

thrust  =  thrust2; 
return ; 

> 
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else 

{ 

if  (  t  <  t2  )  /**  missile  flight  time  **/ 

{ 

w  =  w  -  (mbrm  *  dt); 
thrust  =  thrust2; 
return ; 

) 

else 

{ 

printf(”\n  missile  burn  time  completed” ) ; 
w  =  w; 
thrust  =  0; 

} 

} 

} 


/***  if  user  select  a  long  range  missile  ***/ 


if  (  mselect  ==  3  ) 

{ 

s  =  S3 ; 

if  (  t  ==  0  )  /***  missile  time  of  fire  ***/ 


{ 


w  =  wlr; 

thrust  =  thrust3; 
return ; 


> 

else 

{ 


if  (  t  <  t3  )  /***  missile  time  of  flight  *•**/ 

{ 

w  =  w  -  (mbrl  *  dt); 
thrust  =  thrust3; 


return ; 


) 


else 


printf(”\n  missile  burn  time  completed"); 
w  =  w ; 
thrust  =  0; 
return ; 
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m__guide(  ) 

/*  This  function  computes  the  missile  command  angles,  alpha 

(missile  angle  of  attack)  and  beta  (missile  side  slip  angle), 
during  a  proportional  navigation  flight  path.  The  command 
angles  are  with  respect  to  the  missile  stability  and  are  the 
output  of  this  function.  */ 

{ 

int  c ; 
int  i , j ; 
float  m [3] ; 
float  t  [3]  [3] ; 
float  losdt [3] ; 
float  closv; 
float  x , y , z ; 


x  =  xtar  -  xmsl; 
y  =  ytar  -  ymsl; 
z  =  ztar  -  zmsl; 


rng  =  sqrt(  x*x  +  y*y  +  z*z  ); 
closv  =  (  orng  -  rng  )/  dt ; 


calc  direction  cosine  of  LOS  in  global  coordinate 


v[0]  =  x/rng; 
V [ 1 ]  =  y/rng; 
v[2]  =  z/rng; 


/**  calc  the  rate  of  change  of  LOS  between  missile  and  target  **/ 


for  (  c  =  0;  c  <  3;  ++c  ) 

( 

losdt[c]  =  (  v[c]  -  ov[c]  )/  dt; 
printf("\n  losdt[#d]  =  " , c , losdt [c] ) ; 

} 


/***  matrix  elements  ***/ 


t[0] 

C0]  = 

t  [  0  ] 

[i]  = 

t  [0] 

II 

i — i 

C\l 
i _ 1 

1 - 1 

1 - 1 

-P 

II 

i — i 

ts* 

i _ i 

1 - 1 

1 _ 1 

■P 

[i]  = 

1 - 1 

1 _ 1 

-P 

n 

i — i 

CM 

i _ i 

1 - 1 

CM 

i _ i 

-P 

II 

i — i 

IS* 

1 _ 1 

t  [2] 

II 

1 — 1 

1 — 1 

1 - 1 

C\J 

i _ i 

■P 

II 

1 — 1 
CM 

i _ i 

cos(gamm)  * 
cos(gamm)  * 
sin ( gamm ) ; 
sin( etam ) ; 
-cos( etam ) ; 
0; 

-sin(gamm) 
-sin ( gamm ) 
cos( gamm ) ; 


cos(etam) ; 
sin( etam) ; 


*  cos(etam); 

*  sin(etam); 
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/***  transform  missile  position  from  global  coordinates  to 

missile  stability  coordinates  ***/ 

for  (  i=0;  i<3;  ++i  ) 

{ 

m[i]  =  0; 

for  (  j  =  0;  j  <  3;  +  +  j  ) 

m[i]  =  m[ i ]  +  t [i ] [ j ] *losdt [ j ] ; 

) 

/***  new  direction  cosine  of  LOS  converted  to  old  direction 

cosine  of  LOS  for  the  next  iteration.  ***/ 

ov [0]  =  v  [0] ; 
ov[1 ]  =  v[l ] ; 
ov [2]  =  v [2 ] ; 

/***  computed  range  between  missile  and  target  is  now  the 

old  range  for  the  next  iteration.  ***/ 

orng  =  rng ; 

printf("\n  old  range  =  iof  meters"  ,  orng  )  ; 

/***  calc  the  missile  command  angles,  alpha  and  beta.  ***/ 

alpha  =  ( ( ko*closv*m[2] )*w)  /  ( 0 . 5*cna*s*rho*vmsl *vmsl ) ; 

beta  =  ( ko*closv*m[ 1 ] )^w  /  ( 0 . 5*cna*s*rho*vmsl *vmsl ) ; 

printf("\n  alpha  =  ^f", alpha); 
printf(M\n  beta  =  $f",beta); 

return ; 


} 
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